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ABSTRACT 

A  new  approach  to  the  stochastic  analysis  of  general  com- 
partment models  is  presented.   The  analysis  is  based  on  the  concept 
of  diffusion  approximations.   The  state  of  a  compartment  system  is 
represented  as  the  superposition  of  a  deterministic  process, 
characterized  by  a  system  of  ordinary  differential  equations,  and 
a  random  noise  process  characterized  by  stochastic  differential 
equations.   All  transition  rate  parameters  are  permitted  to  be  time 
dependent.   Numerical  solutions  are  presented  for  the  two-compartment 
case,  and  compared  with  those  of  Cardenas  and  Matis  (1975a) . 
Extensions  to  non-linear  compartment  models  are  discussed. 
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1.   Introduction 

Quantitative  representation  of  the  distribution  over  time 
of  a  drug  or  pollutant  in  a  human  or  animal  body  is  now  often 
carried  out  in  terms  of  compartment  models .   Such  elements  as  the 
blood,  gut,  liver,  lean  tissue,  etc.,  are  characterized  as  homo- 
geneous compartments,   within  which  the  drug  resides  for  a  time, 
later  to  transit  to  another  compartment,  perhaps  recycling  but 
eventually  vanishing  because  of  evacuation  or  metabolism  mechanisms. 

Classical  papers  in  this  area,  e.g.  that  by  Bischoff  et  al . 
(1971)  concerned  with  use  of  methotrexate  in  cancer  therapy,  pro- 
posed deterministic  descriptions  of  flows  between  compartments, 
and  stocks  within,  using  differential  equation  systems.   However, 
variations  in  drug  concentrations  over  replicated  experiments 
suggest  a  probabilistic  or  stochastic  representation.   Stochastic 
compartment  analysis  has  undergone  rapid  development  for  several 
years,  and  might  be  categorized  as  follows.   First,  many  papers 
have  emphasized  the  mathematical  formulation  and  exploration  of  the 
stochastic  behavior  of  a  variety  of  compartment  models;  papers  of 
this  type  are  exemplified  by  Bright  (1973),  Matis  (1972,  1974, 
1975a,  1975b),  Marcus  (1975),  Matis  and  Hartley  (1972),  Purdue 
(1974a,  1974b,  1975) ,  Rescigno  (1973)  ,  Rubinow  and  Winzer  (1971)  , 
Thakur,  Rescigno,  and  Schafer  (1972,  1973),  and  Thron  (1972). 
Second,  several  authors  have  investigated  problems  of  statistical 
inference  for  the  transition  rates  in  stochastic  compartment  models; 
see  Burkinshaw  and  Marshall  (1971) ,  Cornfield,  Steinfeld,  and 


Greenhouse  (1970),  Kodell  and  Matis  (1976),  Matis  and  Hartley  (1972), 
Rodda  et  al .  (1975),  and  Shah  (1976).   Last,  the  actual  application 
of  compartment  models  to  biomedical,  ecological,  and  chemical  systems, 
as  well  as  to  pharmacokinetics,  has  been  carried  out,  for  example, 
by  Metzler  (1971),  Rodenbeck  et  al.  (1975),  Sheppard  (1962),  Siegel, 
Cooper,  and  Meisner  (1968)  and  Uppuluri  and  Bernard  (1967) .   A  large 
bibliography  is  given  in  Jacquez  (1972) . 

The  present  paper  falls  into  the  first  category:   we  suggest 
a  new  approach  to  stochastic  compartment  analysis  based  on  mathe- 
matical diffusion  processes ;  see  Feller  (1971)  for  an  introduction, 
and  Arnold  (1973)  and  Gikhman  and  Skorohod  (1971)  for  a  systematic 
treatment.   In  short,  we  write  stochastic  differential  equations 
to  describe  compartment  concentration  changes,  and  show  that  in  a 
natural  limit  interesting  joint  distributions  are  jointly  Gaussian 
or  Normal.   In  the  literature  a  variety  of  different  compartment 
models  have  been  formulated  for  one  up  to  n-compartment ,  reversible 
or  irreversible,  open  or  closed,  systems.   These  have  generally 
been  characterized  by  discrete-valued  state  variables  for  compart- 
ment contents,  the  latter  changing  according  to  multivariate 

birth-and-death  processes.   Kolmogorov  forward  equations  for  the 
transition  probabilities  are  then  derived,  and  solved  by  transforms. 

Purdue  (1974)  notes  the  similarity  of  these  models  to  those 
for  infinite  server  queues.   We  develop  and  discuss  this  connection 
in  the  appendix  of  this  paper.   For  such  discrete  models  it  is  possible 
to  write  down  expressions  for  the  moments  of  compartment  contents  at 


a  fixed  time;  under  certain  conditions  (Poisson  inputs,  transition 
rates  depend  linearly  on  compartment  contents)  the  contents  of  the 
compartments  are  statistically  independent  and  Poisson  distributed. 
Cardenas  and  Matis  (1975a ,b)  have  extensively  investigated  such 
models.   The  approach  taken  here,  that  of  replacing  a  discrete  birth 
and  death  process  with  a  continuous  diffusion,  is  approximate.   How- 
ever, the  results  agree  with  the  exact  solution  in  the  case  of  linear 
transition  rates  up  to  second  moments.   Furthermore,  the  limiting 
process  may  be  shown  rigorously  to  be  normal,  allowing  the  use  of 
statistical  inference  procedures  validated  for  normal  distributions 
and  processes.   Perhaps  more  important  is  the  fact  that  the  diffusion 
approximation  applies  to  situations  with  nonlinear  transition 
Michaelis-Menten-type  rates,  where  "exact"  discrete-state  methods 
yield  cumbersome  results;  see  McNeil  and  Schach  (1973)  for  various 
examples  concerning  populations.   It  is  reassuring  to  find  that  the 
diffusion  approximations  often  agree  exceptionally  well  with  the 
birth  and  death  solution;  see  Gaver  and  Lehoczky  (19  75,19  76)  for 
numerical  comparisons. 

2.   A  Model 

We  consider  a  general  n-compartment  model  with  time-dependent 
transition  rates.   Let   C± (t) ,  i  = 1, . . . ,n   represent  the  contents 
of  compartments   l,...,n   respectively  at  time   t.   We  assume  that 
the  state  of  the  system  C(t)  =  (^ (t) , . . . ,Cn (t) )   is  a  Markov  process 
with  transition  probabilities  given  by 

P(A±(t)=l,  Aj(t)=0,  j^i|C(t))  =  NXQi(t)dt  +  o(dt) 

i  =  1 ,  .  .  .  ,n 

P(Ai(t)=-l,  A.(t)=Of  j  f    i|C(t))  =  XiQ(t)  C±(t)dt  +  o(dt)      ^^ 

i  =  1 , . . . ,n 
P(A.(t)=-l,  Aj(t)=+1,  Ak(t)=0,  k  +    i,j|C(t))  =  X±j(t)  Ci(t)dt+o(dt) 

for   1  £  if  j  £  n,  i  ^  j#  X^  1  0  • 


where   A. (t)  =  C. (t  +  dt)  -  C. (t) ,  and  N   is  a  parameter  that  can 
be  thought  of  as  a  dose  level,  eventually  allowed  to  be  large. 

The  first  of  the  three  transition  probabilities  represents 
an  increase  in  compartment   i   of  size  1  from  the  outside.   The 
second  represents  a  decrease  in  compartment   i   of  size  1  to  the 
outside.   Finally,  the  third  represents  a  transfer  of  size  1  from 
compartment   i   to  compartment   j . 

The  above  formulation  is  rather  general,  and  allows  many 
different  compartment  systems  to  be  studied  simultaneously.   It 
allows  for  an  open  system  if   Xn.(t)  >  0   for  some   i  =  l,...,n 
or  a  closed  system  if   A...  (t)  =  0   for   i  =  l,m,n.   The  flows  can 
be  either  reversible  if   X . . (t)  >  0   implies   X . . (t)  >  0   or 
irreversible  if   X . . (t)  >  0   implies   X . . (t)  =  0 . 

We  wish  to  study  the  case  where   I-_-i  c-  (t)   is  large,  say 
proportional  to  some  number   N.   A  diffusion  approximation  arises 
when   N  becomes  large  and  we  expand   C(t)   into  terms  of  order  N 
and   /N.   We  can  insure  that   I-_i  C. (t)   is  large  in  two  different 
ways.   For  open  systems  we  assume,  as  given  in  (1) ,  that  the 
transition  probabilities  into  the  system  from  the  outside  are  them- 
selves directly  proportional  to   N.   For  closed  systems  where  there 
is  no  input  from  the  outside,  we  must  assume  that  each  compartment 
has  contents  proportional  to  N   at  time  0,  that  is   C. (0)  =  Nc. (0) . 

The  effect  of  assuming   X-_-i  C. (t)   to  be  proportional  to   N 
is  that  transitions  of  all  types  occur  very  frequently.   As  a  result, 
in  any  short  time  interval  there  will  be  many  transitions  of  each 


type.   Each  type  of  transition  represents  a  change  of   +1   or   -1 
in  the  contents  of  a  certain  compartment.   In  a  short  time  interval 
the  total  change  will  be  a  sum  of  independent  Bernoulli  random 
variables,  and,  as  such,  normally  distributed  by  virtue  of  the  central 
limit  theorem.   Actually,  much  more  is  true.   Focusing  on  the  com- 
partment process   {C(t)  ,  t  >_  0}   over  time  rather  than  at  a  single 
fixed  time  we  see  that  the  process  will  have  normally  distributed 
increments  and  hence  be  Gaussian.   This  observation  is  the  key  to 
understanding  the  diffusion  approximation  approach.   The  mathematical 
development  of  this  approach  is  not  given  here;  see  Kurtz  (19  71) 
and  Barbour  (1972) ,  and  the  paper  by  Schach  and  McNeil  (1973) . 


3.   Mathematical  Development 

In  this  section  we  derive  the  diffusion  approximation  approach 
outlined  above  in  a  manner  that  seems  intuitively  appealing.   The 
exposition  will  be  in  terms  of  (Ito-type)  stochastic  differential 
equations;  although,  as  will  appear  subsequently,  the  stochastic 
differential  equations  describing  stochastic  variations  in  compart- 
ment contents  are  not  ambiguous  of  interpretation.   The  derivation 
to  be  presented  next  is  supplemented  by  further  discussion,  reserved 
for  the  Appendix,  that  relates  the  present  theory  to  the  approach 
by  Barbour  (1974),  and  also  to  infinite  server  queueing-type  models. 


Let,  then,  dCi(t)  =  C^t  +  dt)  -  (^(t)   represent  the  change 
in  the  contents  of  compartment   i   in  time   (t,  t  +  dt) .   The 
change   dC . (t)   may  be  viewed  as  a  sum  of  independent  random 
variables : 

(a)  inputs  from  outside  the  system,  with  mean  and  variance 
NXQi(t)  dt; 

(b)  inputs  from  compartment   j   (j  ^  i) ,  with  mean  and  variance 

A  .  .  (t)  C.  (t)  dt; 
ji      j 

(c)  outputs  to  compartment   k   (k  7*  i)  ,  with  mean  and  variance 

A.,(t)  C. (t)  dt; 
lk      1 

(d)  outputs  from  the  system  via  compartment   i,  with  mean  and 

variance   A.~(t)  C.  dt. 
i0      1 

Represent  the  change  in  contents  of  compartment   i   in  the  manner 
of  diffusion,  i.e.  in  terms  of  drift  and  diffusion  terms: 

n  n 

dC.  (t)  =  [NA..(t)  -   7  A..(t)  C.  (t)  +   J"  A..  C.(t)]  dt 
1  Oi       .**,  in      l       .L.     ]i   i 

3=1  J         :=i  J   J 


n   

+  /NlT-  dWn.  (t)  -   J"   /A..(t)  C.  (t)  dW.  .  (t)         (2) 

Ol      Ol  Ln  11         1  1] 


n 
+  J      /Aj;L(t)  Cj(t)  dWj±(t) 

i  =  1,2, ...,n;  C(0)  =  Nc(0) ,  and   {W..(t),  0    i,  j  <  n,  i  ?   j} 
is  a  family  of  standard  Wiener  processes ,  see  Arnold  (1974)  . 


Expression  (2)  is  motivated  as  follows:   the  drift  term,  in 
brackets,  represents  the  expected  change  in   C. (t)   over  (short) 
time  interval   (t,  t  +  dt) ,  given  the  states  of  all  compartments 
and  the  corresponding  infinitesimal  flow  rates  or  transition 
probabilities.   The  remaining  terms  represent  the  various  random 
components  of  change,  modeled  as  Gaussian  random  variables;  the 
latter  is  plausible  as   N  -»■  *>   on  central  limit  theorem  grounds. 
For  instance,  arrivals  from  outside  the  system  in   (t,  t  +  dt) 

are  Poissonian,  and  hence  for  large   N,  nearly  Gaussian, 
with  mean  and  variance  specified  by  (a)  above,  and  this  accounts  for 
the  first  bracketed  and  unbracketed  terms  in  (2) . 
It  is  shown  in  Kurtz  (1971)  that 

C.(t) 

— - ►  c.  (t)    as   N  ->  oo   in  probability;      (3) 

see  also  Barbour  (1974) .   By  analogy  with  the  central  limit  theorem, 

consider  the  normalized  process 

C.  (t)  -  Nc.  (t) 

x.  (t)  =  — —  ;  (4) 

1  /N 


hence  express  the  concentration  as 


C.  (t)  =  Nc.  (t)  +  )/N  X.  (t)  .  (5) 


In  matrix  form,  C(t)  =  Nc(t)  +  /N  X(t) 


The  vector  function   Nc(t)   is  referred  to  as  the  deterministic 
approximation,  and  /n  X(t)   is  a  stochastic  noise  process  superimposed 
upon  the  deterministic  approximation.   It  remains  to  characterize 
c(t)   and  X(t) .   The  stochastic  process   X(t)  =  (X, (t) , . . . ,  X  (t) ) 
defined  by  (4)  satisfies  a  stochastic  differential  equation  similar 
to  (2),  and  it  can  be  derived  from  (2)  by  an  application  of  Ito's 
Lemma  (Arnold,  p.  90,  or  Gikhman  and  Skohorod,  p.  24).   The  stochastic 
differential  equation  governing   X(t)   is  given  by 


dX(t)    =  A(t)    X(t)    dt   +   B(t)    dW(t) 


-    /N    (c'(t)    -    Xn(t)    -  A(t)    c(t))    dt      (6) 


where 


A(t)    =   a.     (t)     ,  1    <    i,    j    <    n 

AJ  1  J 


n 


with      a.  .  (t)    =    A  .  .  (t)       for      i-^j,a..(t)=-      T      A..(t) 
ID  3i  ?    J         ii  j£0      ID 


j^i 


^0(t)    =    (X01(t) X0n(t))T 


rJ 


w(t)   =   (w01(t) ,...,wQn(t) ,  w10(t) ,...,wn0(t) ,  w12(t) ,... ,wln(t) , 

.. .,   w         , (t) )T 
n,n-l 


B(t)    =    (bik(t))  l<i<n,    1    <   k   <    n(n+l) 


with 


b..(t)  =  /Xn. (t)  i  =  l,...,n 


iiv  '      Oi 


bi(n+i)(t)  =  "  /AiO(t)  Ci(t)/N 


bi(2n+(i-l)  (n-l)+j-l)  (t)    "bj (2n+(i-l)  (n-l)+j-l)  (t) 


=  "  ATj  (t)  Ci(t)/N  ,     j  >  i 


bi(2n+(i-l)  (n-l)+j)  (t)    'bj (2n+(i-l)  (n-l)+j)  (t) 


=  -  /A  .  (t)  Ci(t)/N  ,       j  <  i 


and  all  other   b.,  =0. 

lk 

We  now  let   N  -*"  +  °°   in  (6)  .   First,  the  terms  in   B(t) 


involving  /A . . (t) C. (t)/N  converge  to   /A . . (t)  c. (t) .   Second, 
equation  (6)  contains  a  term  proportional  to   N.   The  coefficient 
of  this  term  must  be  identically  0  for  all   t  >_  0   or  (6)  will 
not  converge  to  a  finite  limit.   Consequently   c(t)   must  satisfy 


c'(t)  =  A  (t)  +  A(t)  c(t)  (7) 


and,  given  that   c(t)   satisfies  (7),  equation  (6)  becomes 


dX(t)  =  A(t)  X(t)  +  B(t)  dW(t)  (8) 

where  C.(t)/N  terms  are  replaced  by  c.(t)  in  ^(t) .  Note  that  the 
modelling  ambiguity  sometimes  associated  with  stochastic  differential 
equation  representations,  see  Arnold  (1973),  Chap.  10,  does  not  arise 
because  the  diffusion  coefficient,  B(-),  of  (8)  does  not  involve   X. 


The  deterministic  approximation  to   C(t)   is  given  by  Nc(t) 
where   c(t)   is  characterized  by  (7).   This  system  of  differential 
equations  is  stable  in  the  sense  of  Liapunov  because  all  eigenvalues 
of   A(t)   have  strictly  negative  real  parts.   The  stochastic  noise 
process  superimposed  on  the  deterministic  approximation  is  given  by 
»^JX(t)   where   X(t)   is  characterized  by  (8).   The  stochastic 
differential  equation  given  by  (8)  describes  a  nonstationary  multi- 
variate Ornstein-Uhlenbeck  process.   Much  is  known  about  such 
processes,  see  Arnold  (1974);  we  will  return  to  a  discussion  of  the 
noise  process  in  a  later  section. 


4 .   The  Deterministic  Approximation 

The  system  of  differential  equations  given  in  (7)  along 
with  an  initial  configuration  c(0)  can  be  easily  solved.  If 
^,0(t)  =  0/  the  closed  system  case,  then 

t 
c(t)  =  exp(/   A(s)ds) -c(O)  (9) 

2 
where   exp(M)  =  I  +M+M/2!  +  •••   for  any  square  matrix  M. 

For  the  open  system,  nonhomogeneous,  case  with  ^0(t)  ^  0   one 

t 
first  multiplies  by   exp(-  /   A(s)ds)   to  obtain 

^  (exp(-/   A(s)ds) -c(t))  =  exp(-/   A(s)ds) • XQ (t)  .  (10) 
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This  equation  can  be  integrated  and  the  result  multiplied  by 

t 
exp(/      A(s)ds)       to   give 
0     ~ 

t  t  t 

c(t)    =    /      exp(f      A(s)ds) *An(u)du   +   exp(/      A(s)ds)-c(0)  (11) 

~  0  u~  ~°  0     "  ^ 


Equation  (11)  provides  a  complete  solution  to  the  problem  of 

describing   c(t)  .   Simplifications  are  not  possible  in  general; 

however,  equation  (11)  can  always  be  solved  numerically  by  standard 

numerical  integration  techniques.   We  can  consider  a  special  case 

which  allows  equation  (11)  to  be  substantially  simplified.   Suppose 

we  assume  that  all  transition  probabilities,  A. . (t) ,  1  <  i  <  n, 

0£j<^n,  i^j,  exhibit  the  same  time  dependence,  or  that 

A. . (t)  =  A. .f(t)   where   f(t)   is  some  function  of  time.   The 
ID        ij 


ma 


trix  A(t)   now  has  the  form 


A(t)  =  Af (t)  (12) 


where 


A  =  (a,..) 


la.   =   A  . 

l  13    :i 

(  a    =  -  I      A.  . 


J  *   i 


j-0   ^ 

Assume  that   A  can  be  diagonalized  and  has  eigenvalues 

6,, 0 «,..., 9  ,  left  eigenvectors   &',...,£  ,  and  right  eigenvectors 
l  2.  n  ^l     wn 

r,  ,...,r  .   This  gives 

~1     ~n        3 
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A(t)  =  RD,(t)L     with    LR  =  I  (13) 


where   R=  (r^,  . .  .  ,rn)  ,  £  =  (^ ,  .  .  .  ^J  \  and   D(t)   is  a 
diagonal  matrix  wii 
integrated  to  give 


diagonal  matrix  with  ith  entry   0.f(t).   Equation  (13)  can  be 


t 
/   A(s)ds  =  RF(u,t)L  (14) 

u 

t 
with   EJ(u,t)   a  diagonal  matrix  having  elements   (0.  /   f(s)ds) 

u 
Finally  (14)  can  be  exponentiated  to  give 


exp  /   A(s)ds  =  RG(u,t)L  (15) 


u 


where   Cl(u,t)   is  a  diagonal  matrix  with  the  ith  diagonal  element 

t 
exp(0.  /   f(s)ds).   Equation  (15)  can  be  plugged  back  into  (11) 

u 
to  give 

t 
c(t)  =  R,  /   £(11,10.1^  <u)du  +  RG(0,t)L-c(0)     (16) 


We  illustrate  the  use  of  (16)  by  considering  an  example 
discussed  by  Cardenas  and  Matis  (1975a) .   In  that  paper  the  authors 
consider  a  case  with  time-dependent  transitions  for  a  closed  two- 
compartment  system.   For  closed  systems  ^0(u)  =  0   and  only  the 
second  term  of  (16)  need  be  calculated.   In  this  case  it  is 
assumed  that 
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and 


X01(t)  =  A02(t)  "  °  ' 

A12(t)  =  .3/(t+l) , 

A1Q(t)  =  .5/(t+l) , 

A21(t)  =  .4/(t+l) 

A2Q(t)  =  .6/(t+l)  . 


A(s)  = 


-.8     .4 
.3     -1 


s+1 


and 


/  -.8     .4 
A  =  ,     f (S)  = 


.3     -1 


s+1  * 


The  matrix  A  has  eigenvalues   8,  =  -.9  +  /. 13 , 

t  e. 

99  =  -.9  -  /.13.   Also   exp(6.  /   f(s)ds)  =  (1+t)   .   It  is  easy 

1  0 
to  calculate 


.12/(.12+a2)   .4a/(.12+a2) 

£(t)  =l  J(  e2  II  2  2  )"JS( 

1  a/4   b/4/\     0     (1+t)  z     \   .12/(.12+bz)   .  4b/  (  .  12+bz)  ' 


a  =  -.1  +  /7T3 


(17) 


b  =  -.1  -  /Til 


Equation  (17)  may  be  expanded  to  yield 
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(18) 


c1(t)  =  (.6386750491c1(0)  +  . 5547001962c2 (0) )  (1+t)  X 

92 
+  (.361324909^(0)  -  .  5547001962c2  (0)  )  (1+t) 

91 
c2(t)  =  (.4160251472^(0)  +  .  3613249509c2  (0)  )  (1+t)  L 

62 
+  (-.4160251472^(0)  +  . 6386750491c2 (0) )  (1+t)  Z 

with 

6,  =  -  .5394448725 

62  =  -1.2605551275 


Transition  probabilities  for  this  model  were  approximated 
by  Cardenas  and  Matis  using  a  truncated  infinite  series  approach. 
These  transition  probabilities  can  be  derived  directly  from  (18) 
by  establishing  appropriate  initial  conditions.   For  example,  by 
setting   c, (0)  =  1,  c~(0)  =  0,  c, (t)   and   c~(t)   become   p, , (t) 
and  p12(t),  respectively.   By  setting   c, (0)  =  0,  c.iO)  =  1, 
c,(t)   and   c2(t)   become   p„, (t)   and   p22(t),  respectively.   The 
results  derived  from  (18)  are  compared  with  the  approximations 
of  Cardenas  and  Matis  in  Table  1.   These  results  illustrate  the 
accuracy  of  the  deterministic  approximations  in  (18)  ,  since  the 
two  answers  agree  to  at  least  8  significant  figures.   Furthermore, 
the  numbers  derived  from  (18)  always  exceed  the  Cardenas  and 
Matis  results.   This  is  reasonable  since  the  latter  approximation 
is  derived  by  truncation  and  is  hence  an  underestimate  of  the  true 
value. 
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The  results  in   (18)   illustrate  another  very  important 
point.   Compartment  modelling  has  been  criticized  because  it  only 
allows  levels  which  are  mixtures  of  exponentials  in   t  when, 
in  fact,  these  models  often  do  not  fit  data  well   (See  Marcus 
(1975)  p.  339  and  Matis  (1972) .    Much  more  complicated   functions 
are  required  to  give  adequate  fits  including  gamma,  Pareto, 
Pareto-exponential,  and  first  passage  density  functions.   One  can 
see  that  all  of  these  types  of  functions  can  be  produced  by 
appropriate  choice  of   f(t).   Mixtures  of  exponentials  arise  when 
f (t)  =  1;  however,  an  enormous  class  of  level  functions  can  be 
produced  by  varying  choices  of   f (t) .   The  introduction  of  time 
dependent  transitions  not  only  adds  realism  but  produces  models 
which  can  fit  real  data  sets. 
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t  Equation  (18)  Cardenas  and  Matis 

4 -term  approximation 

1  .5902421829  .5902421826 

2  .4435620915  .4435620872 
P13_(t)   3  .3652902915  .3652902729 

4  .3155676392  .3155675836 

5  .2807030964  .2807029846 

1  .5151767471  .5151767467 

2  .3596617846  .3596617805 
p22(t)   3  .2823115397  .2823115195 

4  .2356325933  .2356325396 

5  .2041834275  .2041833201 

1  .1501308716  .1501308715 

2  .1678006137  .1678006133 
P2i  (t)   3  .1659575075  .1659575068 

4  .1598700918  .1598700879 

5  .1530393378  .1530393292 

1  .1125981537  .1125981539 

2  .1258504603  .1258504600 
p12(t)   3  .1244681306  .1244681301 

4  .1199025689  .1199025659 

5  .1147795033  .1147794969 

TABLE  1:  A  Comparison  of  the  Deterministic  Approximation 
(18)  with  the  Cardenas  and  Matis  Four  Term 
Approximation . 
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5 .   The  Noise  Process 

The  stochastic  elements  of  the  compartment  process  (C(t),  t  >_  0} 
are  embodied  in  the  nonstationary  multivariate  Ornstein-Uhlenbeck 
process  described  by  (8) .   We  summarize  a  few  key  results  about 
such  processes.   These  results  are  well  known  and  available  in 
Arnold,  Section  8.2.   First,  we  assume   X(0)  =  0   that  is  the  initial 
condition  is  nonstochastic  and  embodied  in   c(0).   The  marginal  dis- 
tribution  of   X(t)   is 


X(t)  ~  */?  (0,  Z(t))  (19) 

*j  *n  /*»  *v 


that  is  an  n-dimensional  multivariate  normal  distribution  where   £ (t) 
is  the  unique  nonnegative  definite  solution  of  the  matrix  Riccati 
equation 

£' (t)  =  A(t)  I (t)  +  E(t)  AT(t)  +  B(t)  BT(t)         (20) 


Since   C(t)  =  Nc(t)  +  /N  X(t)   we  find 


C(t)  ro   77  (Nc(t)  ,  NE(t))  (21) 


Let   B(t)  BT(t)  =  D(t)   and   D(t)  =  (d..(t))   with 


n  n 

1  x   (t)  c  (t)  +  : 

j=0   1]      1       y 


d..(t)  =  (An.  (t)  +  I      A   (t)  c.(t)  +  I      A   (t)  c.(t)) 
11        ui       .Q      i]     1      i=1   j  j-     J 


(22) 


d±j(t)  =  -  (A±j(t)  Ci(t)  +  Aj±(t)  Cj(t)) 
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The  covariance  matrix   E (t)   can  be  calculated  by  various 
numerical  methods  applied  to  (20)  .   We  illustrate  the  calculation 
in  the  case   n  =  2   studied  by  Cardenas  and  Matis. 


A(t)  = 


-U10(t)+A12(t))  A21(t) 


x12(t)       -(x20(t)+x21(t)> 


£(t)  = 


Equation  (20)  becomes 


axl(t)    a12(t) 


a12(t)    a22(t) 


where 


s(t)  =  G(t)  s(t)  +  H(t)  (23) 


s(t)  =  (a11(t)/  c12(t)/  a22(t))' 


G(t) 

A/ 


-2(X10(t)  +  X12(t))  2A21(t)  ° 

X12(t)       -(  X10(t)  +  X20(t)  +  A12(t)  +  X21(t))         X21(t) 
0  2X12(t)  ~2(  A20(t)  +  X21( 

XQ1(t)  +  (X1Q(t)  +  X12(t)c1(t)  +  X21(t)c2(t)) 

-  (  X12(t)c1(t)  +  X21(t)c2(t)) 
X02(t)  +  x12(t)cl(t)  +  U20(t)  +  X21(t)c2(t)) 


Equation  (23)  can  be  solved  in  exactly  the  same  manner  as  (7)  with 

initial  conditions   s(0)  =  0.  that  is   C(0)   is  deterministic.   We 

find 
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t     t 

s(t)  =  /   exp(/   G(s)ds) -H(u)du  .  (24) 

mt  0  U   ~  ** 


For  the  case  considered  by  Cardenas  and  Matis 


+  kS 


The  matrix   G   has  eigenvalues   <j>.  ,  0O,  <|>_   where 

ts*  1     2.  i 

<k    =    26^^    =   -1.978889745 

d>      =      9      +9      =-18. 
y2  12  ' 


/ 

-1.6 

.8 

0 

get)  - 

i+t  I 

.  3 

-1.8 
.6 

-2 

^3  =  2e2  =  -2-521110255 


e  -i        e  -i 

kl;L(u+l)  x   +  k12(u+l) 

V1     .     ,    ^V1 


H(u)  =1  k21(u+l)      +  k22(u+l) 


22 

e,-i         e  -l 

k31(u+l)  X   +  k32(u+l) 


where  the   k. .   coefficients  are  computing  from  (18).   The  other 
factor  in  the  integrand  is  given  by 


0 


/tg(s)ds  =  R     0       (££) 


2 


exp  I      G(s)ds  =  K^     u        llT^        u       l 


with 
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1  1 

.6513878188      -.25 
.4243060905      -.75 


.407957184  .7085463502 
.4615384615  -.3076923077 
.1305557202     -.4998540425 


s(t)   can  now  be  easily  found.   Each  of  the  three  variance-covariance 


terms  is  a  linear  combination  of  terms  of  the  form 


(J) .       (J)  .  —  (j) . 
(1+t)  ""((l+t)  :   x  -  l)/(c})j-(J)i)  . 


The  Ornstein  Uhlenbeck  process  described  by  (8)  also  has  a 
known  covariance  function   k(s/t)  =  Cov(X(s),  X(t)).   The  function 
k(s,t)   is  given  by  the  equation   (T   denotes  transpose) 

min(s,t)    -,  ,     _     _, 

k(s,t)  =  £(s)     /      *  x(u)  B(u)  B/(u)  ($  -"(u))1  du  $x(t)   (25) 

with 

u 
$ (u)  =  exp(/   A(s)ds)  . 
0 

We  do  not  use  the  covariance  function  in  this  paper;  however, 
this  function  is  important  in  the  development  of  a  statistical 
analysis  of  data  from  a  compartment  model.   One  may  have  data,  say 
readings  of  drug  concentrations  in  the  blood  stream,  and  wish  to 
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estimate  the  parameters  of  the  model,  the   X. .'s.   Given  the  data 
have  a  multivariate  normal  distribution,  the  likelihood  function 
can  be  written  and  estimators  derived.   The  papers  of  Hartley 
and  Matis  (1971)  and  Kodell  and  Matis  (1976)  provide  analyses 
in  special  cases.   The  results  of  this  paper  indicate  that  these 
kinds  of  results  can  be  carried  much  further  into  far  more  general 
compartment  models. 


6 .   Steady  State 

In  the  case  ^(t)  7*  0,  an  open  system,  the  system  will 
approach  steady  state  if   X.  .  (t)  ->  X .  .   as   t  -*■  +  ».   Steady  state 
is  a  condition  whereby  the  probability  distribution  which  describes 
the  marginal  distribution  of  the  system  becomes  time  homogeneous 
even  though  the  actual  contents  of  the  compartments  continue  to 
vary  according  to  (8) .   The  steady-state  solution  for  a  simple  example 
is  described  in  the  Appendix. 

The  process   {C(t),  t  >  0}   in  steady  state  can  also  be 
approximated  by   Nc  +  /N  X(t)   where 


0  =  Xn  +  Ac    or    c  =  -A  1Xn  (26) 

and   (X(t),  t  >  0}   satisfies 
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dX(t)  =  AX(t)  +  B  dW(t)  (27) 

a  stationary  multivariate  Ornstein-Uhlenbeck  process.  This  means 
that 

£(t)^9?n(-NA"1X  ,  NE)  (28) 


where   Z   is  the  unique  nonnegative  definition  solution  of 

AI  +  ZAT  =  -BBT  .  (29) 

As  an  illustration  we  consider  a  two  compartment  open  system, 
either  reversible  or  irreversible 


"(X10+X12)         X21 


A  = 


o     -(^o+^i* 

X01+(A10+A12)C1+X21C2 
S8*  ~  I       "  <A12C1+  X21C25   !  (30) 

X02+X12Cl+(X20+X21)c2 


X20+X21     X12   \/X01 


c    = 


X21  X1Q  +  X12/\X02 

(X20+X21} (X10+X12)_X12X21 
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)/2(A10+A12)  2X21  0      -1 

=  I     X12      (X10+A20+A12  +A21>      "X21    I  S5"    <31) 
°22      ^     °  -^12        2(X20+X21> 


Equations  (30)  and  (31)  thus  complete  the  description  of 
the  multivariate  normal  density  given  in  (28) .   While  the   n  =  2 
case  can  be  done  in  closed  form,  the  general  case  can  be  easily 
solved  by  numerical  methods. 


7.   Extension  to  Nonlinear  Models 

The  technique  of  diffusion  approximations  outlined  and 
applied  to  general  linear  compartment  models  in  the  previous 
section  also  provides  a  tractable  approach  to  handling  nonlinear 
models.   The  compartment  modelling  literature  is  vast  and  the 
area  is  still  in  rapid  development.   Nevertheless,  papers  involving 
a  stochastic  process  approach  assume  a  linear  system,  that  is, 
one  in  which  transition  rates  from   i   to   j   are  proportional  to 
the  contents  of  the  ith  compartment.   It  is  unlikely  that  real 
pharmacokinetic  processes  operate  so  simply.   Rather,  the 
transition  rates  are  more  likely  to  be  complicated  functions  of 
the  contents  of  several  compartment  levels,  and  of  time.   Such 
complicated  models  have  not  appeared  in  the  literature,  perhaps 
because  they  present  serious  mathematical  intractabilities  using 
the  standard  Kolmogorov  forward  equations  analysis. 
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The  method  of  diffusion  approximations  can  effectively 
handle  nonlinear  transition  rates.   The  deterministic  approximation 
will  be  a  system  of  nonlinear  ordinary  differential  equations. 
This  will  only  rarely  be  solvable  in  closed  form.   Nevertheless, 
the  numerical  solution  of  such  systems  is  a  standard  topic  in 
numerical  analysis.   The  important  observation  is  that  the  stochastic 
differential  equation  representing  the  noise  process  superimposed 
upon  the  deterministic  process  will  be  nonstationary  but  linear 
(often  Ornstein-Uhlenbeck) .   This  means  that  the  diffusion  approxi- 
mation approach  will  yield  analytic  results  for  nonlinear  systems, 
since  many  results  are  readily  available  for  linear  stochastic 
differential  equations  (see  Arnold,  Chapter  8) .   It  is  hoped  that 
this  approach  will  encourage  pharmacokinetic  researchers  to  con- 
sider introducing  nonlinear  models  should  an  analysis  of  data 
indicate  the  need. 

Two  additional  points  should  be  made.   First,  the  analysis 
presented  can  be  easily  carried  out  when  the  process  is  represented 
in  terms  of  volume  and  concentration  rather  than  in  terms  of 
total  contents.   Second,  the  method  of  diffusion  approximations 
illustrates  that  the  compartment  system  contents  will  have  a 
marginal  Gaussian  distribution.   As  such,  one  is  in  a  position 
to  give  a  statistical  analysis  of  such  a  system  either  to  estimate 
transition  rates  or,  indeed,  to  design  an  experiment  to  make 
proper  inferences  about  such  parameters. 
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APPENDIX 

The  purpose  of  this  appendix  is  to  further  justify  the 
earlier  analysis,  in  particular,  equations  (7)  and  (8)  and  their 
implications.   We  do  this  in  certain  tractable  simple  cases,  namely, 
those  of  one  and  two  components. 

1.   The  Theory  of  Kurtz  (1971)  and  of  Barbour  (1974) . 

In  Section  2  of  Barbour  (1974) ,  hereafter  called  B, 
assumptions  are  given  that  imply  convergence  of  the  sequence  of 
processes 

/C  (t)        \    C  (t)  -  Nc(t) 
yM(t)  =  /N  [^-—   -  c(t)   =  ^ ~ f  N  =  l/2f3 (A_1} 


^N  \   N      'O 


/N 


to  the  diffusion   y(t),  as   N  -*■  °°   in   D(0,T),  i.e.  the  space  of  all 

to 

right-continuous  functions  with  left-hand  limits.   Note  that  the 

y(t)   of   B  will  turn  out  to  be  the  same  as  our  stochastic  noise 
a 

process   X(t).   We  shall  verify  these  assumptions  for  illustrative 
cases. 

First  consider  the  sequence  of  continuous  parameter  Markov 
processes   xN(t)  =  C  (t)/N   on   [0,T],  where   C N(t)   is  the  many- 
compartment  model  whose  transitions  are  described  by  our  (1) . 
Notice  that  the  changes  in   xrT(t)   are  of  size  0,  +  1/N,  and  that 
the  infinitesimal  transition  rates  are  of  the  form  NA^.,  and 
NX. .(C. (t)/N),  as  specified  by  (2.1,B).   We  now  discuss  the  assump- 
tion and  the  Theorem  K  of  B  for  these  illustrative  situations. 
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Single-Compartment  Model.   Let   An, (t)  =  An, ,  and 
A10(t)  =  A10   constant  and  positive,  and   n  =  1,  which  represents 
a  single-compartment  system. 

It  is  easily  seen  that  there  is  a  unique  function   £(t) 
satisfying,  for   0  <_  t   <_   T,  the  natural  deterministic  equation, 
(2.2,B) : 

V   (t)=  AQ1  -  A1Q£(t)  •  (A"2 


The  solution  is,  putting   £(0)  =  c(0), 

#m     ~XWt         A01     M        "^lO1,          /      ,n,          X01    \       "^lO*          A01  ,A    I 

£(t)    =  c(0)e  +  t (1-e  )    =      c(0)    -  t e  +  r ,    (A-3 

A10  \        A10  /  A10 

which  is  unique;  hence  B's  Assumption  A  is  justified.  Suppose 
0  <  c(0)  <  A  /A  for  example;  other  cases  may  be  treated  in 
analogous  fashion.   Then,  in  the  terminology  of  Barbour  (1974) , 

S  =  |c(0),  c(0)  e"XT  +  ^i  (1  -  e"Xl°T)J  (A-4) 


and 


S£ 


[A        —AT 
c(0)  -  e,  c(0)  e"XT  +  ^  (1-e   10  )  +  el  ,      (A-5) 


where   0  <  e    <    c(0)   suffices.   Clearly  the  transition  rates  are 
multinomial  within   S£,  and  B's  Assumption  B  is  also  justified.  Now 
Theorem  K  follows,  and  implies  that   YN(t)  =  v^N  (x^(t)  -  £(t)) 
converges  to  the  diffusion   y(t),  y(0)  =  0,  with  the  characteristic 
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iey(t) 

function   $(6,t)  =  E[e       ]   satisfying  the  differential  equation 


If  =  -6Aio  If  "  I  e2[Aoi  +  X105(t)1*        (A"6) 


recognizable  as  describing  the  Ornstein-Uhlenbeck  process 
whose  stochastic  differential  equation  is 


dy  =  -A1Qydt  +  /XQ1  +  A10£(t)  dW  ,  (A-7) 


this  being  the  same  s.d.e.  as  that  obeyed  by  our  noise  process 
X   for  the  present  setup.   Hence  the  procedure  leading  to  our 
equations  (7)  and  (8)  yields  the  same  results  as  those  rigorously 
established  by  Kurtz,  adapted  by  Barbour. 

Note,  too,  that  a  steady-state  solution  will  exist:  simply 
set   £'  =  0   in  (A-2)  to  discover  the  deterministic  component: 
S(°°)  =  Xoi^A10*   The  stocnastic  noise  has,  from  (A-6)  or  (A-7), 
variance  equal  to  the  mean,  namely,  A01/A1().   Thus  the  steady-state 
compartment  content  is,  according  to  our  theory,  approximately 
/^(NA01/X10,  NA01/^10) /  as   N  becomes  large.   This  is  in  accord 
with  the  fact  that  the  exact  distribution  is  Poisson   (NA   /A   )  , 
as  is  well-known,  and  shown  again  subsequently  in  this  Appendix. 
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Two-Component  Model.   Here   C(t)  =  (C^t),  C2(t)),  the 
allowed  transitions  and  rates  are  summarized  below. 

Transitions  Transition  Rates 

cx(t) ,  c2(t)  + 

a)  Cx(t)  +  1,  C2(t)  NXQ1 

(c1(t)\ 

b)  c1(t)  -  i,  c2(t)  nx10  y—x—) 

c)  C1(t)  -  1,  C2(t)  +1  NX 


(A-8) 

(Cl{t)) 

12  \   N   / 


d)  C1(t)  ,     C2(t)  -  1  NX 


20 


State  scaling,  as  in  (2.1,B)  alters  the  +  1  changes  to  +  1/N. 
The  function   £(t)  =    (£,  (t)  ,  £0(t))   satisfies  the  differential 


equations 


q  -  x01  -  (x10  +  x12)?1 


(A-9) 
"12s!         ,v20s2 


£A    -    Xn  0£..     -    Xor.£ 


For  simplicity,  we  have  not  allowed  entries  into  compartment  2 
from  outside,  and  have  also  omitted  back  flow  from  compartment  2 
to  compartment  1.   The  two  equations  are  readily  solved  to  produce 
a  unique  solution  for  all   t  >_  0 : 

q(t)  =  c1(0)  exp[-U10  +  X12)t]  +  j —^ [1  -  exp(-(X10+  X12)t)J 
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52(t) 


C20  e^P(-X20t)  + 


01       12 


X10+  A12 


20 


[1  -  exp(-X2Qt) 


0,(0)  -  x 


01 


10+  X12 


12 


A20  "  (X10  +  A12} 


x     [exp(-(X10  +  X12)t)  -  exp(-A2Qt)]        (A-10) 


Thus  Assumption  A  is  satisfied,  with  c.(0),  c2^°)  >  °*   Assumption  B 
is  also  satisfied  with 


where 


0  <  e  <  minlj^CT),  £2  (T)  ]  (A-ll) 


£.  (T)  =     min     £  (T)  >  0,    i  =  1,2; 
0  <  t  <  T 


It  is  clear  from  the  differential  equations  that  such  values  will 
always  exist.   Hence 


fc(t,..^  (^-c,t.) 


converges  weakly  to  a  bivariate  diffusion.   Proceeding  further, 
evaluate  the  characteristic  function  of  the  limiting  noise, 

E[exp(i91y1(t)  +  i62y2(t))]  =  $ (B± , 02 , t) ;      (A-12) 

The  latter  satisfies  the  following  partial  differential  equation, 
from  (2.7,B) : 
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•  ■  =   [-e^o  +  i-Q1  +  e2)x12]$0i  -  e2*20%2 

-   \    [e2lX01  +  9lX10^1(t)  +  (-ei+e2)2\2^(t)+  62X20^2(t)]$ 

(A-13) 

This  is  recognizable  as  describing  a  bivariate  Ornstein-Uhlenbeck 
process.   Now  one  may  differentiate  (A-13)  at   6,  =   62  =  0, 
utilizing  the  fact  that 

Vi(t)  =  E[y?(t)]  =  -*Q  e  (0,0,t),     i  =  1,2, 


and 


v12(t)  =  E[Yl(t)  y2(t)]  =  -$e  Q  (0,0, t) 


to  derive  differential  equations  for  the  elements  of  the  variance- 
covariance  matrix  of   y.   Not  surprisingly,  these  differential 
equations  agree  with  (23)  (with   A. .   constant,  and,  in  this  case, 
X21  =  0),  i.e.,  V1(t)  =  axl(t),  V12(t)  =  a12(t),  V2  =  a22(t). 
Thus  if  the  same  initial  conditions  are  adopted  for   £ ,  c   and 
for   y,  and   X   the  deterministic  and  noise  components  of  the 
approximations  of  this  paper  and  of   B   are  seen  to  be  identical. 
The  validation  can  be  extended  with  no  difficulty  in  principle 
to  the  general  situation. 
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2 .   Relation  to  Infinite-Server  and  Linear  Markov  Population  Models 

The  class  of  compartment  models  described  by  the  transition 
scheme  (1)  are  similar  to  the  classical  Poisson-arrival  infinite- 
server  models  of  queueing  theory  and  identical  to  certain  Markov 
population  models  of  Bartlett  (1949 )  ,  and  Kingman  (1969 )  •   We 
now  briefly  explore  the  connection,  illustrating  by  the  single- 
compartment  and  two -compartment  examples.   In  particular,  it  will 
be  shown  that,  even  in  the  time-dependent  parameter  case,  limiting 
Gaussian  marginal  distributions  occur  for  compartment  contents, 
in  agreement  with  the  development  of  the  paper. 

Single-Compartment  Model.  Let   C,(t,u)   denote  the  number 
of  arrivals  to  the  system  that  have  occurred  in  the  time  interval 
(u,t),  0  <_   u  <_   t,  and  that  are  present  in  compartment  1  at  time  t. 
Let   I, (t,u)   be  the  probability  that  an  arrival  occurring  at 
time   u   is  present  in  compartment  1  at  time   t.   Consider  the 
generating  function  of   C,(t,u), 

C,  (t,u) 
E[z11      ]  =  g(z1,t;u)  (A-14) 

Now  by  the  assumption  of  Poissonian  arrivals,  see  (1) ,  we  may 
write 

g(z1,t;u) =  [z1I1(t,u)NX01(u)du+  1-  I1 (t ,u) NXQ1 (u) du]  g (z± , t ;u + du) , 

(A-15) 
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where  the  bracketted  expression  is  the  generating  function  of  the 
contribution  to   C,(t,u)   from  arrivals  during   (u,  u  +  du)  ; 
independence  permits  the  multiplcation .   Subtract   g(z, ,t;u  +  du) 
from  both  sides,  divide  by   du,  and  let  du  •*•  0.   The  differential 
equation 

dg(z, ,t;u) 
5u =  V^i't''^     [(^-D  Ix(t,u)  NAQ1(u)dul      (A-16) 

results,  the  solution  of  which  is,  when   u   is  set  equal  to  zero, 

C,  (t)  t 

g(z1,t;0)=  E[z1-L    |C1(0)=0]  =  exp  [  (z^l)  N  /  I±  (t  ,u)  \Q1  (u)  du]  , 

(A-17) 

the  generating  function  of  the  Poisson  distribution  with  mean 

t 

N  /   I1(t,u)  XQ1(u)du.   It  now  follows  immediately  from  the  central 
limit  theorem  that  the  distribution  of 


t 
C1(t)  -  N  /   I1(t,u)  AQ1(u)du 

X1(t)    = (A-18) 


converges  weakly  to  the  normal  law  with  mean  zero  and  variance 

t 

/   I1(t,u)  AQ1(u)du.   Furthermore,  it  is  easy  to  verify  that 


Var[y(t)]  =  /   i  (t,u)  A_.du  =  -°i  (i  -  e   10  ) 
0    x       U1      A10 
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when   AQ1   and  \^Q    are  both  constants,  and   y(t)   is  the  solution 
of  the  stochastic  differential  equation  (A-7)   which  results  from 
following  the  approximation  program  that  is  the  topic  of  the  paper; 
we  start  with   (^(0)  =  0.   Verification  in  the  time-dependent 
parameter  case  is  also  possible,  but  is  more  difficult.   Thus  the 
correct  marginal  distribution  is  verifable  directly  for  this, 
the  simplest,  case.   With  more  labor  the  joint  distribution  at 
different  times  may  be  found  using  an  extension  of  the  technique. 


Two -Compartment  Model.   Consider  the  generating  function 
of   (C,(t,u),  C2(t,u)),  denoted  by   g (z± , z2 , t ;u) .   Then  an  argument 
analogous  to  that  joint  presented  shows  that   g   satisfies  a 
differential  equation,  which  when  solved  yields 


g(z1,z2,t;0)  =  exp[N{(z1-l)  /   I1(t,u)  XQ1(u)du 


t 

0 


t 
+  (z9-l)  /   I9(t,u)  X   (u)du)]      (A-19) 
0 

implying  that  the  contents  of  the  two  compartments  are  Poisson 
distributed  at  time   t,  and  that   C^t)   and   C2(t)   are  independent 
The  same  central  limit  theorem  argument  now  shows  joint  normality 
as   N  ■>  °°   in  agreement  with  the  claims  of  the  body  of  the  paper. 
Calculations  omitted  here  show  that  the  diffusion  approximation  pro- 
duces first  and  second  moments  that  agree  with  those  of  the  Poisson 
input  model  above.   In  particular,  the  correlation  term,  a-,-,  of  (20) 
equals  zero,  signifying  the  independence  that  we  noted  in  (A-19) . 
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